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In developing a flexible body spacecraft simulation for the Launch Abort System of the 
Orion vehicle, when a rapid mass depletion takes place, the dynamics problem with time 
varying eigenmodes had to be addressed. Three different techniques were implemented, 
with different trade-offs made between performance and fidelity. A number of technical 
issues had to be solved in the process. This paper covers the background of the variable 
mass flexibility problem, the three approaches to simulating it, and the technical issues 
that were solved in formulating and implementing them. 


Nomenclature 

a/ Inertial linear acceleration of the origin of the body reference frame 

a n Inertial linear acceleration of node n 

a c n Inertial acceleration of center of mass of nodal body n 

(20 Inertial linear acceleration of the mass center of the body 

C Damping constant matrix of the body 

c n Position vector of the mass center of nodal body n with respect to node n 

F Total force on the body 

F ex t , n External force on nodal body n 

I n Moment of inertia of nodal body n about node n 

70 Moment of inertia of nodal body n about its center of mass 

K ee Elastic stiffness of body 

in Number of mode shapes used to describe flexible deformation 
m n Mass of nodal body n 
rriT Total mass of the body 

M ee Elastic mass matrix of the body, symmetric (m x m) 

M re Rigid-elastic coupling matrix of the body (6 x to) 

M rr Rigid mass matrix of the body, symmetric (6 x 6) 

N Number of nodal bodies 

q Array of modal coordinates of the body (m x 1) 

r n Position vector of node n with respect to body reference frame after deformation 
T n Total torque on nodal body n about node n 
Text,n External torque on nodal body n about node n 
70 Total torque on the body about its mass center 

u n Displacement of node n due to deformation with respect to body reference frame 
x f Inertial position vector of the origin of the body reference frame 
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x n Inertial position vector of node n 

x 1 ^ Inertial position vector of center of mass of nodal body n 

a f Inertial rotation angles of body reference frame for small rotation 

<p n ,r Shape function of node n for its r th mode of deformation 
'i/’n.r Slope function of node n for its r th mode of deformation 
T Modal transformation matrix 

£ Transformed modal coordinates 

3 Reference frame of the body 

Reference frame of nodal body n, located at node n 
p n Position vector of node n with respect to body reference frame before deformation 
ujb Inertial angular velocity of the body when treated as a rigid body 

u f Inertial angular velocity of the body reference frame 

u) n Inertial angular velocity of nodal body n 

V Matrix expression for a vector V in cross product with another vector 
[] T Superscript T represents transpose of a matrix 


I. Introduction 

In support of NASA’s Constellation Program, 1 a series of pad abort flight tests had been planned to 
demonstrate the ability of the Launch Abort System (LAS) to carry the Crew Exploration Vehicle (CEV), 
or Orion, and its crew to safety in case of a mishap during launch. The first of these flight tests, known as 
Pad Abort 1 (PA-1), 2 has taken place, and it verified the basic functionality of the LAS from the pad based 
on the preliminary design (see Fig. 1 for the PA-1 test vehicle configuration). 


Pad Abort 1 Test Configuration 



Figure 1. Pad Abort 1 configuration of the CEV+LAS system. 


To support analysis, it had been decided to increase the fidelity of the original rigid body PA-1 Guidance, 
Navigation, and Controls (GN&C) simulation, developed at the NASA Lyndon B. Johnson Space Center 
(JSC), by adding structural flexibility to the simulation to account for the deflections of the abort and 
attitude control motor jets as well as the inertial sensors. In some scenarios, these effects are important 
factors in determining ascent abort trajectories. 

The equations of motion (EOMs) for a flexible body are traditionally formulated in terms of a set of 
assumed modes and six rigid degrees of freedom (DOFs). Free vibration mode shapes are the natural choice 
for the assumed modes in the fixed mass case, but for variable mass, vibration mode shapes change as the 
mass of the system depletes. We are not aware of there being a standard way of treating this variable mass 
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flexibility problem, with the literature on the subject quite sparse. 

From the implementation standpoint, the most straightforward approach to the variable mass flexibility 
problem is to continue using a set of fixed assumed mode shapes for integrating the EOMs. We call it 
Method 1, and it is essentially the standard Assumed-modes method, 3 except that the mass matrix and 
other parameters in the EOMs are now time-dependent. This approach was used by Tobbe et al. 4 

Given that the eigenmodes change with mass depletion, we had sought to characterize the accuracy of 
using a set of fixed assumed modes as the basis for integrating the EOMs. This characterization was pursued 
on two fronts: (1) To develop stand-alone tests of how accurately one reproduces the evolving eigenmodes 
by solving the eigenvalue problem in this fixed basis as the mass depletes; and (2) To develop two alternative 
techniques where the modes used in integrating the EOMs accurately reflect the evolving eigenmodes at all 
times. We then compared the latter two simulation approaches with the first one as to the performance and 
the expected degree of fidelity. In the end, this allowed us to select the higher performance approach as our 
analysis tool, after showing that it agreed well with the slower, but a priori higher fidelity ones. We call the 
two additional approaches Method 2 and Method 3. 

Using continuously varying eigenmodes would introduce eigenmode derivative terms into the EOMs which 
would be very cumbersome to compute numerically. It is, however, sufficient to use any fixed set of mode 
shapes as the integration basis, provided that it spans the same subspace as the chosen number of eigenmodes. 
To increase fidelity, the integration basis can be updated periodically to reflect the evolution of the subspace 
spanned by the chosen number of time-dependent eigenmodes. 

For the PA-1 problem, system level Nastran data for the combined CEV/LAS system is available for 
the fully fueled state, fully expended or dry vehicle state, and two intermediate mass states. If we only 
switched between the four available sets of modes during the PA-1 flight trajectory, re-expressing the modal 
state vectors in the updated integration basis could introduce discontinuities in the physical state vectors at 
transitions. Instead, in our Methods 2 and 3 we generate accurate updated eigenmodes during the run by 
solving the vibration eigenvalue problem in a larger basis than the number of modes used in integrating the 
EOMs. This can be done as frequently as needed, lessening the possible discontinuities. Methods 2 and 3 
are distinguished by the type of the larger fixed basis they use. 

The paper is organized as follows. Section II is devoted to the flexible body equations of motion. Sec- 
tion III presents the three methods of simulating the variable mass problem, the numerical results comparing 
the three methods, and the necessary background. This section is deliberately qualitative. Section IV con- 
sists of three subsections devoted to the solution of three technical issues raised earlier in the paper. It also 
contains numerical results on mode accuracy. Section V summarizes the conclusions of the paper. 

II. Flexible Body Dynamics Equations of Motion 

There is a large amount of literature on both dynamics of rigid bodies losing mass and dynamics of flexible 
bodies of constant mass. However, flexible body effects of systems losing mass, such as rockets, have been 
treated in detail by only a few authors. A derivation of equations of a rigid rocket with internal flow that can 
be extended directly to include flexibility of the rocket is provided in Meirovitch. 5 Meirovitch 6 also derives 
the partial differential equations of a flexible rocket with internal flow and presents a closed form solution of 
a special case. Banerjee,' which lists a few other works on dynamics of flexible bodies losing mass, provides 
a new formulation for a variable mass flexible body system, which is suitable for finite element-based model 
of a rocket. Tobbe et al. 4 also address variable mass flexible structures in support of the propellant depletion 
for the Ares I launch vehicle, another element of the NASA Constellation Program. 

We have used a formulation based on the central ideas of Meirovitch 5 and Banerjee 7 in deriving the 
EOMs of the combined CEV/LAS system for PA-1. The effect of internal flow of gases on rigid and flexible 
motion has not been taken into account since it is considered to be small compared to the magnitude of 
the propulsion forces. Also, no data is available to account for this effect, and therefore inclusion is beyond 
the scope of this work. With these approximations, the effect of the ejecting masses becomes equivalent to 
jet thrusts applied at the jet nozzles for the current configuration of the system. Moreover, the so-called 
geometric stiffness considered in Banerjee' has not been included, because the system is considered to be 
sufficiently stiff for the motion of the rocket. 

The motion of the flexible body is described by the six degree of freedom (DOF) motion of a reference 
frame 3 coincident with the structural reference frame of the body, and a small flexural motion of the body 
with respect to this frame. Craig 3 expresses the flexural part of the motion as a linear combination of a set 
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Figure 2. Nodal body description of a flexible body. 


of assumed vibration modes of the system. A Finite Element Model (FEM) of the CEV/LAS system was 
obtained from NASA’s Constellation program database; it considers the system as a collection of rigid nodal 
bodies. Figure 2 depicts the flexible body schematically. 

The linear displacement u„ and rotation 9 n of node n with respect to the reference frame is expressed in 
terms of a set of shape and slope functions 0 n ,r and ij) n}r and the modal variables q r as 

m m 

u n = <Pn,r(lr and 0 n = '^2'lp n , r q r (1) 

r=l r= 1 

Since the flexural motion is considered to be small, the small angular motion 9 n is treated as a vector. 
In Equation (1) it is assumed that for the purpose of our analysis flexural displacements can be adequately 
expressed in terms of m shape and slope functions. Selection of the shape and slope functions is discussed 
in the next section. 


II. A. Kinematics 

The linear velocity and acceleration of node n with respect to frame are obtained by differentiating u n 

m m 

U n =^2<Pn,rqr and U„ = ^ <Pn,rQr (2) 

r= 1 r= 1 

Similarly, the angular velocity and acceleration of node n in frame Q is obtained by differentiating 9 n 

m m 

0 n = ^ ^ ^n,rQ.r Bud 0 n = ^ ^ ^n^rQr (3) 

r= 1 r= 1 

From Fig. 1 we can see that r n , the position of node n in frame ^5 and its inertial position x n are given 

m 

T n = Pn H - = Pn “I - ^ ^ Pn,rQr (4) 


by 


and 


C n = Xf + f n = Xf + p n + (pn, r qr 


(5) 


r—1 


Using kinematics, we get the expressions for inertial linear and angular velocities as well as linear and 
angular accelerations of node n. For inertial acceleration of node n we have 


= a/ + Qf x r n + u n + uif x (ujf x r n ) + 2 ujf x 


( 6 ) 
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Inertial angular velocity and angular acceleration of node n are given by 

io n = LOf + 9 n and io n = u)f + 0 n + ujf x 0 n (7) 

Since nodal body n is rigid, the inertial acceleration of its mass center is given by 

CL n — (l n T Ut n X C n T LOn X (^n X C n ) (8) 


II. B. Dynamics 

The inertia force F* and inertia torque T* of the nodal body n about its mass center are given by 


= — 777 n c 

x n " l "n Uj n 


(9) 


and 

T* = — ( I n U) n + U n I n LO n ) (10) 

For each nodal body, the sum of the inertia force, forces that acts on it from other nodal bodies, and 
any external force that act directly on it, must be zero. When we sum the inertia and applied forces on all 
the nodal bodies, the internal forces between nodal bodies cancel out according to Newton’s Third Law of 
Motion, resulting in the translational EOM 

E^* + E^ = ° ( n ) 

n l 


where F ext j is the total external force acting on the flexible body at node l. 

Taking the moment of the inertia forces and applied forces on the nodal body n about the origin of the 
reference frame 9 and adding them for all the nodal bodies, and using the fact that moments of the forces 
between the nodal bodies cancel out as a consequence to Newton’s Third Law, we get the rotational EOM 

E^n + 50 X^+tl+^fiX F exU + Y, Text, l = 0 (12) 

n l l 


where T extt i is the external moment acting on the body at node l. 

For the flexible body EOM, we consider the virtual work done by the inertia and external forces and 
moments, and the elastic forces between the nodal bodies due to virtual flexural deformation of the system. 
Expressions for the virtual translation and rotation of application points of inertia and external forces is 
obtained from Eq. (1). The virtual work consideration leads to the flexible body EOM 


E 






m 

■F: + - Y K r,i9i 


i = 1 


m 

^ ' C r ,iXli -(- y ' Wl,r • T’ext,/ T y \ l,r • F ex i^ — 0 
i= 1 l l 


(13) 


for r = 1, 2, ..., m, where K and C are the mx m stiffness and modal damping constant matrices, respectively. 
After substituting the expressions for F* and T*, the EOMs of Eqs. (11), (12) and (13) may be written in 
matrix form as 

M rr M re 
M? e M ee 



where Af 



b r and b e are, respectively, (6 x 1) and (?n x 1) arrays of nonlinear terms, i.e., terms 


from F* and T* in Eqs. (11), (12), and (13) that do not contain the acceleration terms a/ , w/ and q\ and 
f r and f e are (6 x 1) and (?n x 1) arrays, respectively, containing the F ext j and T exti i terms. 

The terms M rr , M re , M ee , b r and b e involve constant integrals, commonly known as modal integrals, 
computed from the shape and slope functions and node mass and moment of inertia (summed over all nodal 


5 of 19 


American Institute of Aeronautics and Astronautics 



bodies). Derivation of these terms is straightforward but very tedious and is omitted here for the sake of 
brevity. Defining g r = f r — b r and g e = f e — b e , the EOMs of the flexible body may be written as 


M rr M re j A f 

Mj e M ee | q 


0 

Kq + Cq 


9r 

9e 


(15) 


Equation (15) is identical in form to the expressions given in both Quiocho et al. 8 and MacLean et al., 9 
with the exception of the single body versus multibody notation. 


III. Numerical Solution 

III. A. Fixed, continuously varying, and piecewise constant bases 

Even for very small systems, for which all DOFs can be kept, one has to choose a set of flex modes to 
separate the rigid and flex DOFs. A convenient choice is to solve the vibration eigenvalue problem, which 
produces 6 rigid and 6 N — 6 flex eigenmodes (where N is the number of 6 DOF rigid nodes in the FEM). 
For large systems, only a small number M of the 6 N — 6 flex modes is retained in the EOMs. A common 
approach is to keep M lowest frequency eigenmodes. 

It is intuitively appealing to generalize this to the variable mass case and consider the integration basis 
of M time-varying eigenmodes. Here integration basis refers to the set of flex modes used in integrating the 
EOMs. While conceptually attractive, this basis of time-varying eigenmodes is not easy to use in practice 
due to the mode derivative terms that would appear in the EOMs, which would require computing these 
derivatives numerically. 

However, the intuitive appeal of the integration basis of time-varying eigenmodes may overstate the 
necessity of such a basis. If all 6 N — 6 flex DOFs were used in integrating the EOMs, any basis of 6 N — 6 
modes, combined with the 6 rigid DOFs in the EOMs, would span the same set of possible deformed states of 
the system in various locations and orientations. A fixed basis of this size would then be rigorously equivalent 
to the basis of 6 N — 6 time- varying eigenmodes. 

Even for a reduced set of M modes, a fixed basis of M modes that spans the same M-dimensional 
subspace at time t as M lowest frequency eigenmodes, will approximately continue to do so for an interval 
of time A t, barring possible mode crossings. At time t + At the fixed basis can be updated to match the 
latest set of M eigenmodes, and the updated basis will remain accurate for the next Af period. We have 
now arrived at the idea of a piecewise constant integration basis, i.e., one that remains constant between 
periodic updates. Its importance is that one integrates the EOMs between times t and t + At just like for any 
fixed integration basis, with no additional mode derivative terms in the EOMs, while allowing it to change 
periodically to track the evolution of the eigenmodes. 

An additional computation required with a piecewise constant integration basis is to re-express the flex 
state of the system in terms of the new basis vectors when the basis is updated. When the basis is complete, 
so that M = 6 N — 6, we can impose continuity conditions on the deformations and their velocities. If the 
flex modes at time t are and the mode amplitudes and their time derivatives are rf^\t) and ? 

(2) 

when the flex modes are updated to (pj , the continuity requires that 

6 6N-6 6JV-6 

r=l f = 1 f=l 

6 6N-6 6JV-6 

r— 1 /= 1 /= 1 

where the rigid modes are included on the left hand side of the equation only. The reason the six rigid mode 
amplitudes are absent from the right hand side is because they are not present in the EOMs, where the 
rigid DOFs are represented in a different way than small deformations. However, we are forced to introduce 
them on the left because the size of the deformation vector is 6 N, so 6 N — 6 flex amplitudes will generally 
not be sufficient to represent it once the basis is changed. The physical interpretation of these seemingly 
superfluous rigid mode amplitudes is that matching the inertial states of all finite element nodes upon the flex 
basis change, also requires an adjustment of the inertial state of the vehicle Structural Reference (SR) frame, 


(16) 

(17) 
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i.e. , of the rigid DOFs. Adjusting the rigid body position, orientation, and their corresponding velocities 
every time the flex basis is updated is possible but not desirable. This is because of the additional small 
errors in the inertial states of the individual nodes that would be introduced by the nonlinearity of the 
3D rotations. We developed a special technique to avoid these SR frame adjustments; it is presented in 
Section IV.C. 

When the basis is incomplete, so that M < 6 N — 6, the equalities of Eqs. (16) and (17) can generally 
not be satisfied. Instead, one seeks to minimize the resulting discontinuity, i.e. the differences between left 
and right hand sides of the two equations, in an average sense. The update algorithm that was used is 
also presented in Section IV.C. While the possibility of discontinuities exists and must be accounted for, we 
know from our prior experience with flexible multibody space robotics simulations that in practice it may 
be negligible. 7 8 We have not observed them in the quantities of physical interest that we have examined for 
the two PA-1 simulations we implemented that use piecewise continuous bases. 

We would like to note that we consider it a possibility that in the limit At — > 0, a piecewise constant 
basis of periodically updated M lowest eigenmodes, coupled with re-expressing the modal variables in the 
new basis at every update, is mathematically equivalent to using M continuous eigenmodes with modal 
derivative terms. However, we have not yet verified this hypothesis. 

Of the three basis types we have described fixed, continuously varying, and piecewise constant - only 
the fixed and piecewise constant bases are of practical use, at least for the PA-1 problem. 

Note that implementing a piecewise constant integration basis requires a means of interpolating M lowest 
frequency eigenvectors during the run. To take full advantage of this type of basis, we want to be able to set 
the basis update frequency to an arbitrary value, which generally requires being able to generate the basis 
vectors accurately at any time during the run. For example, in our PA-1 simulation where we were only 
provided with 4 sets of vibration eigenmodes corresponding to 4 different fuel levels, we chose to update the 
integration basis 20 times during the run. 

The next section describes three different ways of treating modes in a variable mass problem that were im- 
plemented for PA-1: a fixed-basis one, and two utilizing piecewise constant integration bases. The differences 
between the latter two are in how mode interpolation is handled. 

III.B. Three methods of the EOM integration 

The available Nastran FEM data for the PA-1 system were for the fully fueled (Full), fully depleted (Empty) 
and two intermediate (Il-intermediate and 12-intermediate) fuel levels; Il-intermediate had half full Abort 
Motors and nearly full Attitude Control Motors, while 12-intermediate had empty Abort Motors and half 
full Attitude Control Motors. Each data set included the first 100 free vibration mode shapes (6 rigid zero 
frequency and 94 flex finite frequency) and their frequencies for the corresponding mass matrix, which was 
also included in the data set. The stiffness matrix was assumed to be independent of the fuel level. Since the 
purpose of the simulation was for GN&C analysis, 34 flex modes were used for integrating the EOMs given 
by Eq. (15), providing frequency coverage up to about twice the maximum frequency relevant for GN&C. 

The three methods that were implemented for integrating Eq. (15) are: 

Method 1: Fixed integration basis 

Method 2: Piecewise constant integration basis - eigenmode interpolation via “ordinary ” large fixed basis 
Method 3: Piecewise constant integration basis - eigenmode interpolation via “ expanded ” large fixed basis 

III.B. 1. Method 1 

In the first method, a fixed set of mode shapes is used as the integration basis. Equation (15) is integrated 
directly, inverting the mass matrix at every time step to obtain the second derivatives from the generalized 
forces. The eigenvalue problem is never solved explicitly at run-time; it is built into the formalism implicitly. 
In practice, a certain number of lowest frequency eigenmodes from one of the mass states are often chosen 
as the basis for state integration across all of the mass states. 

Various modal integrals in Eq. (15), namely M re , M ee , and those appearing in the nonlinear terms on the 
right hand side were computed with the above fixed basis in the four available mass states and interpolated 
for other mass states that occur during the run based on those four sets of modal integrals. 
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An important advantage of performing integration in a fixed set of dynamical variables is that no dis- 
continuities occur in any physical quantity. It is also relatively computationally efficient since one does not 
perform transformations of various quantities between different bases, nor ever solve an eigenvalue problem. 

A drawback of Method 1 is the lack of controllability of the eigenmode modeling accuracy. A basis 
consisting of M exact eigenmodes from one mass state may not yield M accurate eigenmodes in another, yet 
the latter are implicitly used at the point in the integration when the system passes through the corresponding 
mass state. 

With the eigenmodes for the PA-1 system available for four mass states, there are four choices for 
selecting the integration basis of 34 flex modes. To characterize the accuracy of a candidate integration 
basis, we first solved the 40 x 40 vibration eigenvalue problem (6 rigid + 34 flex DOFs) in each of the three 
mass states other than the mass state the basis was drawn from. Using a simple mathematical expression for 
the “distances” between these computed and the corresponding true eigenmodes, we obtained a quantitative 
measure of accuracy with which the candidate basis reproduces the true eigenmodes in various mass states. 
This distance measure is one of the two mode characterization techniques we have used that are presented 
in Section IV. A. 

Table 1 in Section IV. A shows the distances between 34 true eigenmodes in the Empty fuel state and 
the corresponding approximate eigenmodes obtained with the basis from the 1 1-intermediate fuel state. 
From that table, we see that the first 24 approximate Empty flex modes are accurate, while the top 12 
are inaccurate, some quite drastically so. These results are typical. After examining the four candidate 
integration bases via this technique, we selected the Il-intermediate set as the better all around integration 
basis to approximate the true eigenvectors across the whole mass spectrum from Full to Empty. As we can 
see, even this choice is far from perfect. 

One can, of course, always seek to improve the accuracy of a fixed integration basis by increasing the 
number of modes in it. This affects performance because the mass matrix that has to be inverted at every 
time step is larger, and the size of other modal integrals increases as well. Furthermore, the appearance of 
higher frequencies often requires reducing the size of the integration time step. 

III.B.2. Method 2 

To deal with the shortcomings of Method 1, which used a fixed integration basis, and to check its accuracy 
for PA-1 modeling as an integrated simulation, rather than just via stand-alone mode proximity tests, a 
simulation with a piecewise constant integration basis was developed. We call this Method 2. 

Method 2 uses an auxiliary larger fixed basis of size Ng + 6 (six rigid + Ng flex modes) for eigenmode 
interpolation needed by the piecewise constant integration basis during the integration basis update. When- 
ever the integration basis needs to be updated, one solves the ( Nb + 6) x (Nb + 6) eigenvalue problem with 
the current mass matrix, retaining Nb < Nb lowest-frequency eigenvectors for use in integrating Eq. (15). 
The eigenmodes are recomputed several times during the run, the state of the system is re-expressed in the 
new integration basis at each update, and the integration continues in the new integration basis until the 
next update. 

An implementation detail worth mentioning is that it can be numerically costly to transform some of 
the higher-order modal integrals from the fixed large basis in which they have been pre-computed to the 
new integration basis at each basis update. It can be more efficient to transform the state vector of the 
system from the integration to the larger fixed basis, compute the nonlinear terms in the fixed basis, and 
then transform them back to the integration basis for actual integration. 

The main advantage of using an appropriately chosen larger auxiliary basis, and then using only a subset 
of its eigenmodes as the integration basis, is that the integration is then performed with accurate eigenmodes 
in all the mass states. While the basis of size Nb + 6 will not yield Nb accurate flex eigenmodes in all the 
mass states, it may yield accurate results for the lowest Nb eigenmodes. This approach can thus be seen as 
a kind of interpolation scheme for the lowest Nb eigenmodes across all the relevant mass states. One can 
also think of it in terms of starting with a larger fixed set of Nb assumed modes in the EOMs of Eq. (15) 
and then performing modal decomposition and truncation to Nb retained modes. 

An important drawback is the potential appearance of discontinuities in the physical quantities important 
for the PA-1 problem at each integration basis recomputation. This may be especially pronounced at mode 
crossings, when what was mode number Nb becomes mode number Nb + 1 and is dropped from the basis 
along with the deformation of the body due to that mode’s former non-zero amplitude. This effect would 
be present even if the basis recomputation were performed at each integration time step. With relatively 
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infrequent updates, mode Nr may jump through several modes to Nr + j with j > 1, and it is even 
possible that more than one mode will leave the original set during a single update. Even without such 
mode crossings, the recomputed set of Nr lowest eigenmodes will span a slightly different subspace of the 
IVA-dimensional basis space, resulting in slight discontinuities that would, however, diminish with the higher 
frequency of integration basis update. 

Valid as the above theoretical considerations may be, past experience with flexible robotic manipulator 
dynamics simulations (such as those in Quiocho et al. 8 ) shows that these discontinuities are not necessarily 
noticeable in practice, at least when updates are frequent enough. For the PA-1 simulations under discussion, 
it appears that whatever discontinuities take place are either small or are not in the quantities that affect 
the integrated GN&C performance. 

The auxiliary basis of Method 2 simply uses the 6 rigid plus all 94 finite frequency modes from one of the 
four available data sets. To assess the accuracy of the four choices the four data sets offer for the auxiliary 
basis, the same kind of stand-alone mode analysis described earlier in Section III.B.l (and presented fully 
in Section IV. A) was performed. Table 3 in Section IV. A shows the accuracy of the 34 eigenmodes in the 
Empty fuel state obtained with the 100 vector auxiliary basis of eigenvectors from the Il-intermediate mass 
state. The 34 eigenmodes are shown to be quite accurate, as they are also found to be in the other mass 
states, so the Il-intermediate basis is a suitable choice of the auxiliary basis. 

III. B. 3. Method 3 

One problem with the basis of Method 2 is its large size relative to the number of retained modes, which 
degrades performance. Another is its “unimprovability” with regard to accuracy. There are four choices 
of the 100 vector auxiliary basis with the four available data sets; they provide whatever accuracy for the 
lowest 34 modes across the mass spectrum that they provide, and once the best basis among the four has 
been chosen, no further improvement in accuracy for the 34 modes is possible. 

To address these problems, we developed a technique for generating a different type of auxiliary basis 
starting with the available eigenmodes in all four mass states. The 6 rigid plus 34 lowest finite-frequency 
modes in one of the sets were used as a ‘starter basis’, and the 34 lowest modes from each of the other 
three sets that were not already accurately represented by the starter basis were analyzed for possible 
addition to that basis. Using an approach inspired by the Gram-Schmidt orthogonalization procedure, 9 10 the 
remaining vectors were analyzed during each iteration, and those that were nearly linearly dependent with 
the starter basis and the other vectors that had already been added to the basis during previous iterations, 
were eliminated. One of the remaining vectors was added to the basis and the iteration then repeated. The 
vector elimination threshold e was set to 0.01, and the resulting 51 vector auxiliary basis (6 rigid plus 45 flex 
modes) provided the accuracy of about 0.01 for the lowest 34 modes across all the mass states. We call this 
type of basis an “expanded” basis, to show that the basis vectors in it are drawn from different mass states, 
unlike what what was done in Methods 1 and 2 above. 

After the 51 x 51 eigenvalue problem is solved, the additional 11 flex modes are not guaranteed to be 
physical and are simply discarded, leaving 34 accurate flex modes for integration. The exact procedure and 
the rationale for it are presented in Section IV. B. 

Stand-alone eigenmode analysis (see Tables 4 and 5 in Section IV. B below) shows that this basis indeed 
produces accurate results for 34 lowest eigenfrequencies and eigenvectors across the whole mass range of the 
simulation. It achieves this accuracy while being considerably smaller in size than the 94 vector basis of 
Method 2. Moreover, the accuracy is more controllable since we can adjust the accuracy threshold e. 

III.C. Numerical Results 

The PA-1 FEM model is quite large (over 70,000 nodal bodies). Mass matrices, mode shapes, and frequencies 
were provided for six rigid modes and 94 lowest frequency vibration modes for each of the four fuel levels, 
which were in turn used to generate four sets of modal integrals. During the simulation run, the flexible 
model remained active for the interval of 11 to roughly 32 seconds, from the time the launch is aborted to 
the time when the LAS separates from the crew module. This approach was used because the LAS is the 
primary contributor to the flexibility of the combined CEV /LAS system. 

Figure 3 shows the linear deflection of one of the Abort Motor jet nozzles in the axial direction for four 
different cases. The 94- modes-1 1-Method 1 shows the data obtained using all 94 mode shapes with the fuel 
level Il-intermediate and Method 1, whereas 34modes-Il-Methodl represents the same mode shapes but 
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keeping only 34 modes while also using Method 1. The 34-modes-I2-Method2 is for the same fuel level and 
number of retained modes but instead using Method 2 with a larger auxiliary basis of 100 modes (6 rigid + 
94 flex). Finally, 34modes-X-Method3 shows the results for Method 3 with 34 retained modes obtained with 
an auxiliary “expanded basis” of 51 modes (6 rigid + 45 flex). 



Figure 3. AM Nozzle 1 flex linear displacements. 


Figures 4 through 6 show selected components of differences in inertial positions of the mass center 
obtained by different methods, with respect to selected methods. Only the components having the largest 
differences are presented. Figure 4 shows the differences in Z coordinates of inertial positions of the mass 
center from solutions using Methods 1 and 2 with respect to Method 3 (considered to be the most accurate 
for a given number of retained modes from theoretical considerations). Results for Methods 1 and 2 are 
labeled 34modes-Methodl and 34modes-Method2, respectively. The difference between the original rigid 
solution and Method 3 is labeled Rigid. 

Figure 5 shows the differences between Y components of inertial positions obtained by Method 1 keeping 
34 modes when mode shapes of different fuel levels are used. The difference is computed with respect to the 
solution that uses mode shapes of the Il-intermediate fuel level. 

Figure 6 shows the differences between Z components of inertial positions obtained by Method 1 keeping 
all 94 modes when mode shapes of different fuel levels are used versus the solution obtained by Method 
3 with 34 retained modes. For completeness, the results also includes the data for the original rigid body 
solution. 

The results of Fig. 3 show that the system is quite stiff and therefore the flexible body trajectories 
computed by the different methods are close to each other. Figures 4 through 6 demonstrate this. Figures 3 
and 4 show that Methods 2 and 3 produce almost indistinguishable results. This is expected since both 
methods use a regularly updated piecewise continuous integration basis, and as Tables 3 and 4 in Section IV. A 
below illustrate, both employ accurate eigenmodes for the instantaneous fueling level during each integration 
basis update. 

The differences in flex deflections on one axis seen in Fig. 3 were found to reduce significantly between 
Method 1 with 94 assumed modes, Method 1 with 45 assumed modes, and Method 2 with 45 retained 
modes. However, the differences in deflections seen with 34 modes do not cause significant differences in the 
computed trajectories. Our interpretation is that 34 modes are sufficient, as expected, for GN&C analysis, 
but some of the other response still requires higher frequency modes for better accuracy. It can also be 
noticed that the rigid trajectory is observably different from the flex trajectories, though the differences are 
not large for the overall distance traveled. 
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Figure 4. PA-1 Mass center trajectory differences for different methods using mode shapes of fuel level Il-intermediate 
with 34 retained modes. Differences computed relative to solution by Method 3 with 34 retained modes. 



Figure 5. PA-1 Mass center trajectory differences using standard assumed-modes method and mode shapes for different 
fuel levels, with 34 retained modes. Differences computed relative to results for fuel level Il-intermediate. 
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Figure 6. PA-1 Mass center trajectory differences using standard assumed-modes method and mode shapes for different 
fuel levels, with 94 retained modes. Differences computed relative to solution by Method 3 with 34 retained modes. 


Based on the results, we conclude that Method 1 with 34 retained modes is sufficient for performing 
Monte Carlo type GN&C analysis runs, because it does not require periodic modal decomposition of the 
system (a computationally expensive task) and is therefore the fastest in execution. 


IV. Technical Developments 


IV. A. Characterization of Mode Shapes 

IV.A.l. Mode proximity measure 


In Methods 2 and 3 we solve the eigenvalue problem in various mass states in specified bases. Method 1 does 
not involve explicit numerical solution of the eigenvalue problem, but vibrations are present in the dynamics 
implicitly, with the implicit basis being the set of assumed modes used plus the additional six rigid body 
modes. 

Suppose we have a set of true and approximate eigenmodes, fame and 4> computed , and both sets are 
normalized with respect to the mass matrix M. We want to find some measure of the magnitude of the 
difference between the true and the corresponding approximate eigenmode to characterize the accuracy of 
the approximate solution. 

Without additional criteria, for example, special emphasis on the accuracy of the displacement of certain 
nodes, a generic measure is needed. A natural distance measure to use in this case is the one based on the 
mass matrix M . The corresponding mode distance criterion is 


D = min 


( ^computed T (frtrue ) M {(ft computed i rue ) 


1 1/2 


(18) 


Note that the eigenvectors produced by different programs for solving the eigenvalue problem may differ from 
each other by a minus sign, so an approximate mode close to the exact one may be pointing in the opposite 
direction. This would result in a large difference vector {<j> computed ~~ fame)- To account for this possibility, 
we compute Eq. (18) for both the sum and the difference of the two vectors, and choose the smallest of the 
two numbers as the accuracy measure. Using the normalization condition on fame and 4> C omputedi we obtain 


D = 


2(l I ^ computed^' fame I ) 


1 1/2 


(19) 
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The maximum value this expression can yield is \/2. when the true and computed vectors are perpendicular. 
The minimum value is zero when the computed eigenvector is identical to the true one. 

For eigenvectors with closely spaced eigenfrequencies, the order of the computed frequencies may be 
opposite to that of the true ones. To analyze for this possibility, we compute the matrix of distances between 
different pairs of computed and true eigenvectors as 


D 


mn 


2 (! ^ I ^computed, m M< t> 


true,n 



(20) 


These values can be recorded as a table and then inspected for possible eigenvector interchanges. 


IV. A. 2. Adequacy of reduced subspace measure 

Another useful measure of what using a reduced subspace does is described here. We take a true eigenvector 
of mass matrix M^ 2 \ expand it in terms of a basis (fin ^ orthonormalized with respect to mass matrix M^\ 

<t>?rue= Y C ^ 1)+ Y Cn ^ ( 21 ) 

n<N r e t n>N rs t 


and observe that the first part of the expansion lies in the reduced subspace spanned by N ret retained 
basis vectors. The second term lies outside that subspace and can never be captured by any approximate 
calculation of the modes restricted to the reduced subspace. Note that we include 6 rigid modes among the 
retained ones. 

Numerically, the weight of the out-of-subspace part relative to the in-subspace part is captured by 


R = 


JX 


n>N ret 


n<N ret 


a 1/2 



(22) 


It is helpful to imagine a vector in a 3-dimensional space, and ask how well it can be represented by a vector 
in a reduced 2-dimensional subspace, i.e., a plane. Referring to Fig. 7, we characterize the subspace adequacy 
by the ratio of the length of the out-of-plane component of the vector, to the length of its in-plane component. 
This ratio is just the tangent of the angle between the vector and the plane, which is the geometric analogy 
for our expression above. 



To evaluate R of Eq. (22), we compute the first N ret coefficients c n from 

n<N ret , (23) 

since we have data for the retained modes. For the remaining c n coefficients, consider <fr ^ M^(jA 2 \ which 
is a square of the norm of the vector <p( 2 ' > with respect to the mass matrix MA\ Substituting the expansion 
(/)<■ 2 ) = '}2n c n ( l ) n' > for 4>^ 2 \ and using the orthonormality relations among <j>n\ we obtain 

♦» T *<V*>-£<4= E'J+E'i < 24 > 

n n<N re t n>N re t 

Since we already know coefficients c n for n < N ret , and since we can compute <j)^ T directly, the 
final expression showing how much of the true vibration eigenvector (jA 2 ^ does not lie in the <pn \ ti < N ret 
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subspace is 


( 2 ) T M«<^ 2 ) - Y. r 

(S n <iV ret C n 

A similar measure using matrix M ^ rather than can also be written. 

IV. A. 3. Numerical results for the accuracy of the PA-1 modes 

Table 1 shows the distances (see Eqs. (18) and (19)) between the computed and true modes in the Empty fuel 
state, the latter obtained with a 40 mode basis of 6 rigid and 34 flex modes drawn from the 1 1-intermediate 
mass state. One can see that the first 22 modes are accurate, while the top 12 modes are not accurate, some 
quite drastically so. Table 2 shows the R values for the Empty eigenmodes with respect to this basis. The 
two criteria are seen to complement each other well. 



i<N r , 


1/2 


1/2 


(25) 


Table 1. Distances D nn between 34 true and computed finite frequency Empty modes, the latter obtained with a 40 
mode (6 rigid -j- 34 flex) basis of eigenmodes from the Il-intermediate mass state. The data is listed row by row. 


0.005229 

0.000334 

0.000763 

0.015117 

0.250930 

1.269756 


0.005729 

0.000370 

0.004134 

0.013731 

0.294966 

1.276961 


0.000112 

0.009187 

0.000743 

0.037520 

0.810824 

0.893421 


0.000361 

0.009899 

0.000139 

0.042375 

0.687222 

1.283446 


0.000269 

0.000866 

0.000494 

0.246295 

0.567872 


0.000024 

0.000473 

0.000501 

0.267012 

1.375331 


Table 2. Measure R for 34 true Empty flex eigenmodes, showing the ratio of the out-of-basis part to the in-basis part 
for each eigenvector with respect to the 40 mode (6 rigid 34 flex) basis of eigenmodes from the Il-intermediate mass 
state. The data is listed row by row. 


0.004599 

0.000358 

0.000804 

0.015179 

0.041344 

0.188542 


0.005130 

0.000388 

0.004596 

0.012181 

0.192263 

0.370136 


0.000124 

0.010180 

0.000700 

0.033963 

0.473543 

0.661967 


0.000403 

0.010884 

0.000130 

0.039985 

0.266537 

0.880828 


0.000320 

0.000811 

0.000394 

0.189347 

0.127798 


0.000026 

0.000505 

0.000412 

0.099630 

0.215277 


Since we chose the number of modes to provide twice the frequency coverage needed for GN&C analysis, 
such basis may still yield adequate simulation results if it is only the higher frequency modes that are 
not accurately represented throughout the whole mass range. The distance measure allows us to test a 
prospective basis by indicating which modes computed with it in different mass states are accurate. 

To further validate the results obtained with this basis, we use a larger auxiliary basis. The larger basis 
is also tested, to verify if it provides accurate results for 34 flex modes in all mass states. Table 3 shows 
the distances between 34 true and computed modes in the Empty fuel state with a 100 vector (6 rigid + 94 
flex) basis drawn from the Il-intermediate mass state. We can see that this basis is reasonably accurate to 
simulate 34 modes with confidence. 

IV. B. “Expanded” Basis 

IV.B.l. Expanded basis algorithm 

We have four sets of PA-1 Nastran data for 6 rigid plus 94 lowest frequency flex modes representing four 
different fuel levels. Choosing one of these sets as the basis for solving the eigenvalue problem in all the mass 
states may yield adequate results, but it is not a controllable procedure beyond the ability to choose among 
the four sets. 
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Table 3. Distances D nn between 34 true and computed finite frequency Empty modes, the latter obtained with a 100 
mode (6 rigid -{- 94 flex) basis of eigenmodes from the Il-intermediate mass state. The data is listed row by row. 


0.001150 

0.000082 

0.000126 

0.004856 

0.003087 

0.018635 


0.001188 

0.000088 

0.002015 

0.010771 

0.014496 

0.012768 


0.000032 

0.002059 

0.000381 

0.006734 

0.018830 

0.011903 


0.000092 

0.002057 

0.000049 

0.005845 

0.014210 

0.008483 


0.000133 

0.000170 

0.000071 

0.017134 

0.003948 


0.000006 

0.000056 

0.000070 

0.005925 

0.015211 


If we hypothetically used all 400 vectors from all of the available mass states as one big basis, the resulting 
eigensolution would provide exact eigenmodes and eigenfrequencies in each of the four mass states. This 
hypothetical basis is unnecessarily large because of many nearly vanishing linear combinations among the 
vectors. We eliminate such near linear dependencies by the Gram-Schmidt inspired procedure with finite 
resolution that is described below. 

If we only need to interpolate a smaller number of modes than 100, this heuristic argument still applies 
if we start with a subset of Nr modes from each of the data sets. 

Step 1: Gather Nr lowest eigenvectors, including 6 rigid modes, from each of the four sets into a set of 

(s) 

4 Nr vectors. Choose one of the four sets of Nr eigenvectors as the seed basis. Call these basis vectors e„ 
and the matrix they come from M^ s \ Set the accuracy threshold e. 

Step 2: Reduce the other 3 Nr vectors by subtracting their components parallel to the seed basis vectors, 

Nr 

v^v' = v — £( v .MW e W) e W. (26) 

n—1 


We use el** that are normalized with respect to their mass matrix M^ s \ Vector v initially stands for any of 
the e„ ^ from the other sets s'. 

Step 3 : Evaluate the magnitude of each reduced v via 

D(v) = (v • M (/) v) 1/2 (27) 


Note that we evaluate each reduced vector’s size with respect to it’s own mass matrix M^ s \ not the M ^ 
of the seed basis. 

Step 4: Eliminate vectors with Z?(v) < e from further consideration. 

If all the reduced vectors v have thus been eliminated, go to Step 7 below to record the final basis, which 
is just the seed basis. 

If any reduced vectors v with the magnitude greater than e still remain, add the one with the largest 
magnitude D(v) to the seed basis. Name it f a . 

NOTE: Do not rescale any of the remaining reduced vectors v. Their unsealed size shows how much of 
an error still remains in the eigenvectors they originate from. 

Step 5 : Subtract the components of the last added basis vector f a from all the remaining non-eliminated 
reduced vectors via 


v = v — 


(v • M«f 0 ) 

(f 0 -M(«)f 0 ) 


(28) 


where the projection is defined with the respect to the seed basis mass matrix M^ s >. 

Step 6: Go back to Step 3 and iterate. 

Step 7: Record the basis. We have generated a basis consisting of the seed basis and a set of vectors that 

(s') 

have been transformed from their original values e„ by a series of subtractions. However, we do not save 

(s') 

the transformed versions, but rather the original e„ vectors - transformations are only used as part of the 
selection process. 
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Step 8: Compute the stiffness matrix among the selected basis vectors. The basis vectors are eigenvectors 
from four different mass matrices M^ s \ so each satisfies 


lif e W=u;W 2 MW e W, « = ,4 (29) 

for its own M^ s \ For two vectors from the same set s, the stiffness mass matrix element between them is 
already known to be 

e m ■ Kt 4 S) = 4 S) 2 S mn (30) 

To compute K among basis vectors from different sets s' and s, we use a mathematical identity involving 
the vibration eigenvalue Eq. (29), 


e 


(«' 

n' 


) 


■ Kei s) 



= OJ, 


(s') 


2 



5 • m (s,) 4 s) 


(31) 


Equations (30) and (31) allowed us to compute the stiffness matrix in the chosen basis even though the 
full stiffness matrix K used in Nastran computations was not available to us. The need to rely on Eq. (31) 
to compute the stiffness matrix elements in this case is why we kept the original, rather than transformed, 
basis vectors in Step 7 above. 

Step 9: We have eliminated linear combinations below e, but linear combinations just above e remain. In 
order to have well-conditioned mass and stiffness matrices, the basis is transformed further. Solving the 
vibration eigenvalue problem 

I<e p = to 2 p M^e p (32) 

in the basis produced in Steps 2-8 results in an orthonormal set of eigenvectors with respect to the seed mass 
matrix M^ s > . This constitutes the final expanded basis. 

Gram-Schmidt orthonormalization procedure could have been used here instead, but this way both 
the mass matrix M ^ and the stiffness matrix K are diagonal in the resulting basis, which is even more 
convenient. 

A procedure based entirely on each eigenvector’s own mass matrix M <s \ rather than the one involving 
the seed basis mass matrix M^ s \ can also be written down. 


IV. B. 2. Accuracy of the expanded basis for the PA-1 problem 

Expanded basis includes flex modes drawn from different mass states. A flex mode that is orthogonal to the 
rigid modes with respect to its own mass matrix will not be so with respect to another mass matrix. To 
ensure that the flex eigenvectors obtained with an expanded basis are orthogonal to the rigid modes, rigid 
modes must be included in the basis generation procedure, and will thus be included in the final basis. 

Tables 4 and 5 show the results for computing Empty modes in terms of the 51 vector basis obtained 
from the seed basis of 40 (6 rigid + 34 flex) Il-intermediate fuel state modes, with the additional 11 modes 
from other sets added by the above algorithm (accuracy set to e = 0.01). 

Table 4. Distances D nn of 34 computed to true finite frequency Empty modes obtained with a 51 mode “expanded” 
basis generated starting from 40 (6 rigid 34 flex) Il-intermediate fueling level modes. The data is listed row by row. 
Distances exceeding the basis selection accuracy of 0.01 are underlined. Data for two additional modes is presented. 
They are seen to be entirely inaccurate. We generally only rely on 34 flex modes to be physical with this basis. 


0.001449 

0.000125 

0.000281 

0.004469 

0.007818 

0.026726 


0.001027 

0.000177 

0.001740 

0.015309 

0.006579 

0.000098 


0.000058 

0.005595 

0.000260 

0.006778 

0.000020 

0.000082 


0.000114 

0.004223 

0.000044 

0.007191 

0.020413 

0.000095 


0.000098 

0.000265 

0.000071 

0.003384 

0.025037 

0.552487 


0.000007 

0.000091 

0.000071 

0.013952 

0.014953 

1.381411 


The results are clearly quite good. As we can see from Table 4, there are several computed modes whose 
distances to the corresponding true modes exceed the threshold value of e = 0.01 that was used in the basis 
selection algorithm. In fact, the expanded basis algorithm’s accuracy threshold e is a necessary, but not a 
sufficient condition, for each computed eigenvector to be within distance e of the true one. Therefore, we 
can expect some of the eigenvectors to be less accurate than e. However, we see that in practice one still 
obtains excellent results. 
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Table 5. Fractional frequency error (/ computed — f true) / / true of 34 computed finite frequency modes in the Empty fueling 
state obtained with a 51 mode “expanded” basis generated starting from 40 (6 rigid -|- 34 flex) Il-intermediate fueling 
level modes. The data is listed row by row. Data for two additional computed modes is also included. 


1.257e-04 1.077e-04 6.045e-08 4.961e-07 1.158e-07 2.046e-08 

3.655e-07 3.200e-07 2.268e-04 1.600e-04 3.881e-07 1.368e-07 


3.134e-07 1.615e-05 

5.980e-06 5.049e-04 


2.532e-07 9.081e-08 -8.321e-09 1.191e-07 

2.212e-05 4.609e-05 6.757e-06 6.123e-05 


6.176e-06 1.204e-05 -2.597e-ll 1.296e-05 1.084e-05 7.687e-06 

1.568e-05 -7.869e-ll -1.362e-ll -8.330e-12 II 4.604e-04 1.345e-03 


IV. C. Modal Truncation of the Equations of Motion 

Here we describe a technique for re-expressing modal variables in the updated integration basis that avoids 
the need for the SR frame repositioning, which, as explained in Section III. A, would otherwise be needed 
when a piecewise constant integration basis is used. We focus specifically on the practically relevant case 
where modal truncation is also present. 

The EOMs in Eq. (15) contain six rigid and m flex states in some fixed basis. In the present case, m is 
larger then the number of modes we wish to integrate. Methods 2 and 3 remove the high frequency modes 
of the changing system periodically, using the procedure described here. We consider free vibration of the 
system without damping for which the motion consists of small rigid and flex displacements. The EOMs 
reduce to 

\m„ m „1 (x,\ , To ol tx,\ _ fol 


where ' *' = l.M “ d ' % “ d S/ ” “ 0,,s “ d r °“ io ” ot body refer ““ fame w “ h 

respect to inertial frame. This equation may be written more concisely as 


MX + I<X = 0 (34) 

where definitions of A, M and K follow from the above equations. Let ’L be the eigenvector matrix of the 
pair [ M , K\ such that M4' = I (identity matrix) and = (w 2 ) (diagonal matrix). Using the special 

property of K that all elements of its first six rows and columns, corresponding to the rigid body modes, 


are zero, it can be shown that 4' and (w 2 ) may be written as f = 


'IVr 'T r 

0 'U 


and (w 2 ) = 


0 0 


The matrices dr and (Co 2 ) have been partitioned according to rigid and flex components. We now truncate 
'T by retaining columns that correspond only up to the desired frequency and define the truncated matrix 
4 / rr . x l , r R 

'T/j = ,r re . Expressing the variables Xf and q by the transformation 

0 ^ eeR 


'LftCft = 


A 'L re R^eR 
*1 'eeR^eR 


Xf 

*1 'eeR^eR. 


We have now determined a transformation for the flex variable q = ^eeR^eR where is the set of new 
modal variables. The subscripts R and e represent the retained and elastic modes, respectively. Note that 
the rigid coordinates are not affected by the transformation. Applying this transformation to Eq. (15) and 
premutiplying the bottom equation by 'T 7 R we get the modally reduced EOMs 


2 C( w fi)Cefl + 


KeRde 


where, M re = M re ^ eeR , M ee = ^ eR M ee ^/ eeR and (u 2 R ) = ^1 eR K^ eeR . In Eq. (36) damping is considered 
to be modal, so that ^ eR C^ eeR = 2 ((tu R ), where C is the modal damping factor. Equation (36) is integrated 
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numerically to propagate the rigid and flex states with time. The system flex variable q is computed from 
£, eR using Eq. (35). 

For a system with changing mass and or stiffness, 'P ee R changes and therefore it requires periodic updates. 
When this happens £ e /j must be re-initialized. The rigid states are propagated without re-initialization. 
During re-initialization it is desired that q ^ where the superscripts W and represent states 

before and after the transition. Using the transformation in Equation (35), it is therefore desired that 
= '^leR^eR- Without modal truncation 'P would be invertible and could be determined 
directly. However, with truncation the dimension of q is greater than that of £ e R and therefore, we can 
only estimate such that it would minimize a norm of the estimation error. One reasonable norm is 

J = i^eeR^eR ~ 9 {1) } T m { e^eR — gd)] where, R is a positive definite matrix. Performing the minimization 
( 2 ) 

with respect to Qr we obtain 

t(2) ,t,(2) T D,T/0) ,t,(2) t p (1) 

S e R ~ eeR ^ eeR\ ^ eeR (37) 

A"( 2 ) is a good choice for R for the estimation of since it would minimize the strain energy corre- 

(2) 

sponding to the estimation error. With that choice, we get the estimate for Qr 

= (^ (2) " 2 )^iS^ (2) 9 (1) (38) 

For reinitialization of we follow the same procedure. In this case M iV is a good choice for the 
positive definite weighting matrix because the norm of the estimation error would be the corresponding 
kinetic energy. We then get 

£r = _1 (39) 

We note that this approach avoids the SR frame repositioning by starting with the eigenvectors 

and zeroing out their rigid parts ’P re _R. Since these transformed vectors are not themselves eigenvectors, M ee 
of Equation (36) is not diagonal in the resulting set of assumed modes. 

V. Conclusions 

We started with the Assumed-modes approach for simulating variable mass flexible spacecraft, using a 
fixed set of mode shapes, but time-dependent mass matrix and other parameters in the EOMs. Given that 
the eigenmodes change with time in such a system, we sought the means of characterizing the accuracy of 
this approach in a systematic way. 

We developed two measures of mode accuracy and applied them to the basis of fixed assumed modes in 
various mass states, allowing us to characterize the ability of such a basis to accurately model the eigenvectors 
at various levels of mass depletion. We found that the Assumed-modes method does not reproduce all of 
the individual evolving eigenmodes accurately at all times, especially the higher frequency ones, making it 
desirable to further assess the quality of the overall simulation results with this method. 

Based on our prior experience with space robotic manipulator simulations, we also developed a simulation 
technique with modal decomposition and truncation in which the set of modes for integrating the EOMs 
is regularly updated to accurately match the evolving vibration eigenvectors. Accurate eigenmodes were 
generated by solving a vibration eigenvalue problem in a larger basis than the number of modes used in the 
EOMs. Two distinct types of such bases were used, resulting in two approaches with modal decomposition 
and truncation. 

These two methods, which agreed well with each other, were shown to reproduce the eigenmodes used in 
the EOMs accurately throughout the whole mass range, thus serving as a benchmark for the Assumed-modes 
method. The comparison showed that the Assumed-modes approach gave good overall results, validating it 
as the high performance analysis tool for the PA-1 GN&C scenario. 

We believe that these three techniques, with different trade-offs between fidelity and performance, offer 
a useful variety of options for analysis as well as for verification and validation purposes. We also hope that 
some of the specific technical problems addressed and solved during the course of this work may be of more 
general interest to the modeling and simulation community. 
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